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I INTRODUCTION 


This  study  describes  and  evaluates  a tracking  algorithm  that  can 
follow  a maneuvering  target  by  filtering  time-sequenced  measurements 
from  multiple  observers.  The  algorithm  can  be  used  when  several  units 
are  passing  targeting  data  to  a single  unit  on  which  the  tracking 
solution  is  computed.  The  data  may  consist  of  (1)  x-and-y  (latitude- 
and-longitude),  (2)  range-and-bearing,  (3)  bearing-only,  or  (4)  course- 
and-speed  measurements,  along  with  estimates  of  the  measurement  errors. 

The  tracking  algorithm  is  applicable  to  the  OTH  (over-the-horizon) 
targeting  situation  shown  in  Figure  1.  Here  two  ships  cross-fix  the 
target  ship  by  using  towed-array  bearing-only  measurements,  and  locate 
each  other  by  range-and-bearing  measurements  using  an  intermediary 
helicopter  as  a radar  target.  The  helicopter  also  relays  the  targeting 
information  from  the  second  ship  to  the  first  ship.  The  tracking 
algorithm  was  specifically  designed  for  this  kind  of  problem,  in  which 
the  locations  of  several  units  (both  friendly  and  target  units)  are 
unknown  but  can  be  deduced  by  time-processing  position  and  velocity 
measurements,  especially  bearing-only  measurements  from  an  ESM  receiver 
or  a passive  sonar.  The  tracking  algorithm  is  called  MURLOC,  an  acronym 
for  Multiunit  Relative  Localization. 

MURLOC  has  several  unique  characteristics ' that  make  it  a versatile 
tracker  in  OTH  targeting  scenarios: 

• All  units  are  correlated  by  employing  a single  state  vector 
that  contains  the  positions  and  velocities  of  all  the  units, 
instead  of  multiple  state  vectors,  one  for  each  unit.  The 
reason  for  using  a super-state  vector  is  to  capture  the 
statistical  correlations  between  units.  Thus  a range-and- 
bearing  measurement  on  an  own-force  unit  automatically 
adjusts  the  estimates  of  relative  position  on  all  units 
including  the  target. 
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FIGURE  1 OTH  APPLICATION  OF  MURLOC 


• Measurements  are  transformed  to  Cartesian  coordinates  before 
they  are  filtered;  thus  MURLOC  uses  a standard  Kalman  filter 
without  having  to  compute  partial  derivatives  at  predicted 
positions,  as  done  in  an  "extended"  Kalman  filter.  This 
"transformed-measurement"  approach  avoids  problems  that  are 
Inherent  in  extended  Kalman  filters. 

• Bearing-only  measurements  are  treated  as  range-and-bearing 
measurements  by  estimating  a pseudo  range  using  a special 
"ellipse  tangent"  algorithm.  The  range-and-bearing  values 
are  then  transformed  to  Cartesian  coordinates  and  filtered. 

The  "transformed-measurement"  approach  to  filtering  and  the 
"ellipse  tangent"  algorithm  for  bearing-only  measurements  seem 
to  work  very  well  in  the  multiunit-type  scenario. 

• A model  noise  covariance  matrix  is  adapted  to  the  measurement 
residuals  so  that  unknown  maneuvers  can  be  tracked.  MURLOC s 
maneuvering-target  adaptive  algorithm  is  simple  and  seems  to 
be  effective. 
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The  impetus  for  creating  MURLOC  was  the  ship-launched  Harpoon 
targeting  problem  of  (1)  processing  data  from  radar,  sonar,  and  ESM 
sensors  that  are  located  on  the  Harpoon  ship  and  other  ships  or  air- 
craft, (2)  tracking  several  targets  and  own-force  units  at  the  same 
time,  and  (3)  producing  decision  aids  to  help  determine  when  to  launch 
a missile.  The  main  computations  needed  are  the  fire-control  solution, 
the  estimate  of  errors  associated  with  the  solution,  the  probability 
that  the  Harpoon  missile  seeker  will  illuminate  the  target,  and  the 
probability  that  the  target  is  within  counterattack  range. 

These  shipboard  computational  requirements  can  be  marginally 
satisfied  with  a multiplicity  of  plotting  techniques,  tables,  graphs, 
slide  rules,  and  hand-held  calculators.  Many  different  techniques  must 
be  developed  for  the  many  different  situations  that  can  arise.  These 
techniques  will  be  complicated  if  they  are  designed  for  situations  in 
which  the  data  are  distributed  in  time  and  come  from  sources  other  than 
own-ship  sensors.  The  calculation  of  targeting  errors  and  other  decision 
aids  will  be  quite  difficult  without  automated  methods. 

An  alternative  to  the  multiplicity  of  techniques  is  to  use  a 
tabletop  programmable  calculator  (such  as  the  Wang  2200  or  HP  9830)  to 
compute  the  required  decision-aid  and  fire-control  parameters.  The 
idea  is  that  all  the  various  targeting  methods  can  be  integrated  into 
one  system  by  developing  programmable-calculator  software  that  can 
incorporate  each  new  measurement  as  it  arrives  in  the  command  center. 

The  data  may  be  from  own-ship  radar,  sonar,  or  ESM;  another  ship's 
sensors;  or  the  LAMPS  helicopter's  sensors.  In  addition  to  originating 
from  different  sources,  the  measurements  may  be  staggered  in  time,  may 
come  in  late,  or  may  be  out  of  time  order.  The  targeting  software 
should  be  capable  of  calculating  tracks  on  several  targets  at  the  same 
time;  and  since  the  position  of  friendly  units  is  important  to  targeting 
solutions,  the  software  should  also  calculate  tracks  on  friendly  units. 

MURLOC  was  created  as  a prototype  algorithm  for  the  tracking 
portion  of  such  a targeting  system.  So  far,  only  the  multiunit,  bearing- 
only,  and  maneuvering-target  tracking  capabilities  of  MURLOC  have  been 
investigated.  A major  unanswered  question  is  whether  or  not  the 
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processing  of  MURLOC's  super-state  vector  requires  so  much  computation 
time  on  a tabletop  programmable  calculator  that  the  tracking  solution 
cannot  keep  up  with  the  data.  The  major  factors  in  this  question  are 
the  number  of  units  in  the  tracking  problem  and  the  data  rate.  Another 
question  is  how  to  filter  data  that  arrive  out  of  time  sequence. 

Because  MURLOC  is  a recursive  algorithm,  solutions  would  probably  have 
to  be  saved  periodically  so  that  the  tracker  could  be  initialized  at 
the  closest  solution  to  the  old  (bui.  newly  arrived)  measurement,  and 
then  recycled  through  the  stored  measurements  following  that  old 
measurement . 

Whether  or  not  MURLOC  is  used  in  a shipboard  targeting  system,  it 
can  stand  alone  as  a method  for  analyzing  multiunit  tracking  scenarios. 
MURLOC  can  be  used  to  process  a time  sequence  of  simulated  random 
measurements  and  predict  the  target  track.  Many  replications  of  the 
random-measurement  sequences  can  be  processed  and  the  average  estimated 
track  compared  to  the  true  track.  This  was,  in  fact,  the  method  used 
to  investigate  MURLOC's  tracking  capability. 

A time  sequence  of  average  position  and  velocity  errors  was 
computed  using  MURLOC  on  50  replications  of  a maneuvering-target  scenario. 
MURLOC  was  compared  to  another  tracking  algorithm  that  used  a batch- 
processed  least-squares  method.  The  comparison  analysis  showed  that 
MURLOC  tracked  the  maneuver  quite  well,  especially  in  estimating  the 
target’s  relative  position.  MURLOC  does,  however,  underestimate  its 
own  errors. 

The  least-squares  algoritnm  is  also  a product  of  the  study. 

Research  on  how  to  perform  weighted  nonlinear  least-squares  computations 
and  how  to  select  the  best  set  of  measurements  resulted  in  a tracking 
methodology  that,  with  further  research,  might  be  useful  in  targeting 
software.  In  any  case,  the  use  of  the  least-squares  algorithm  provided 
a good  benchmark  for  judging  MURLOC's  tracking  capability. 

The  MURLOC  computer  program  is  coded  in  FORTRAN  for  SRI's  CDC  6400 
computer.  The  source  deck  is  approximately  500  cards,  and  with  array 
dimensions  that  can  accommodate  10  units,  the  program  takes  about 
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15,100  words  ot  memory.  The  running  time  depends  on  a number  of  factors; 
for  example,  one  replication  of  a scenario  with  3 units,  10  time  steps, 
and  3 measurements  each  time  step,  required  about  16  seconds  of  computer 
time.  The  MURLOC  computer  program  source  code  is  listed  in  Appendix  A. 


II  MURLOC  TRACKING  ALGORITHM 


Thia  chapter  describes  the  major  features  of  the  MURLOC  algorithm 
and  gives  an  example  of  how  multiunit  measurements  can  be  processed 
with  MURLOC. 


MURLOC  D«»scription 


MURLOC  can  accept  time-sequenced  position  or  velocity  data  on 
multiple  targets  from  multiple  observers.  MUEILOC  processes  the  input 
data  at  each  step  in  time  and  estimates  the  position  and  velocity 
vectors,  the  range  and  bearing  between  units,  and  error  ellipses. 

MURLOC  uses  a Kalman  filter  with  one  multiunit  state  vector  that  causes 
full  correlation  between  units.  A unique  adaptive  covariance  scheme 
reduces  the  problem  of  the  linear  motion  assumption.  Table  1 defines 
the  various  symbols  that  are  used  in  this  chapter,  and  Table  2 summarizes 

th^e  major  equations  that  are  the  essence  of  the  MURLOC  algorithm. 

( 

f 

' 1 . State  Vector 

MURLOC  uses  one  multiunit  state  vector.  For  example,  if  there 
are  five  units  in  the  problem,  the  state  vector  will  have  20  components; 
the  first  four  components  are  the  position  and  velocity  (x  y i y)  of  the 
first  unit,  the  next  four  components  are  the  position  and  velocity  of 
the  second  unit,  and  so  on.  The  covariance  matrix  of  this  example  is  a 
20-by-20  matrix;  thus,  the  state  values  can  become  correlated,  not  only 
by  position  and  velocity  correlation  on  a single  unit,  but  also  by 
correlation  between  units. 

MURLOC  is  designed  so  that  units  report  data  on  other  units; 
no  distinction  is  made  between  friendly  and  enemy  units  other  than 
knowing  which  is  which  by  unit  number.  Thus,  data  from  a radar  fix  on 
a companion  ship  are  processed  in  the  same  way  as  data  from  an  ESM  fix 
on  an  enemy  ship.  This  design  was  chosen  because  of  the  real-world 
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Table  1 


FILTER  VARIABLES 


Variable 

De  f inition 

* 

Dimension 

X 

o 

Estimate  of  the  state  vector,  given  the 

measurement  at  t (the  filtered  state) 

0 

4n-by-l 

p 

0 

Covariance  matrix  of  x 

0 

4n-by-4n 

y 

Measurement  vector  at  time  tj^ 

2 -by-1 

R 

Covariance  matrix  of  y 

2 -by-2 

X 

Estimate  of  the  state  vector  at  tj^  before 
the  measurement  at  tj^  is  filtered  (the 
predicted  state) 

4n-by-l 

p 

Covariance  matrix  of  x 

4n-by-4n 

F 

State  transition  matrix  from  time  t 
to  tj^  ° 

4n-by-4n 

Q 

Model  noise  covariance  matrix 

4n-by-4n 

M 

Measurement  matrix  (state-to-measurement 
transformation) 

2 -by-4n 

r 

Predicted  measurement  residual  error 
(also  called  "Innovation") 

2 -by-1 

s 

Covariance  matrix  of  r 

2 -by-2 

K 

Filter  gain  matrix 

4n-by-2 

^1 

Estimate  of  the  state  vector,  given  the 
measurement  at  t]^  (the  filtered  state) 

4n-by-l 

Covariance  matrix  of  Xj^ 

4n-by-4n 

n = Number  of  units  being  tracked. 
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Table  2 


MURLOC  ADAPTIVE  ALGORITHM 
Filtered  State  at  Time  tg 

X 

o 

p 

o 

Measurement  at  Time  tj^ 

y 

R 

Predicted  State  at  Time  t]^ 

X = F X 

o 

r = y - M X 

S=R  + MFP  f'^m'^ 
o 

P = 1 - exp(-^  r"^  S ^ r) 

T = t,  - t 

1 o 

T 

q = ( Ti  r^  Tj/T  ) 

Q = P q 

p = F p F^  + Q 
o 

Filtered  State  at  Time  tj^ 

K = P M^  (R  + M P M^)‘^ 

L = I - KM 

Xj^  = X + K r 

T T 

P^=LPL  + KRK 


problem  of  tracking  friendly  units,  in  addition  to  the  more  obvious 
problem  of  tracking;  enemy  units.  Because  the  position-velocity  state 
of  friendly  units  is  part  of  the  problem,  the  friendly  units'  state 
vectors  must  also  be  estimated.  The  only  way  this  can  be  done  without 
restrictions  and  ad  hoc  assumptions  is  to  incorporate  all  stale  vectors 
into  one  lar^e  state  vector.  Althouj;h  this  approach  requires  considerable 
computation  to  process  one  measurement,  it  is  beneficial  because  it  uses 
data  properly.  Kor  example,  the  ran^e -and-bear Ing  data  from  Unit  i on 
Unit  5 will,  in  general,  change  the  state  estimates  of  all  the  units 
because  of  past  Interaction  among  them. 

The  state  vector  at  lime  t is  denoted  x . The  estimate  of 

o o 

X is  denoted  x , and  the  covariance  matrix  of  the  error  in  x Is 
o o o 

denoted  P . The  state  at  time  t = t,  evolves  fran  the  state  at  time  t 
o I o 

by; 


X = K X + w 
o 


where  F is  the  transition  matrix  from  t to  t,,  and  w is  a random 
vector  that  is  normally  distributed  with  zero  mean  and  covariance  Q.  As 
discussed  later,  the  Q-iiuitrix  is  defined  as  a function  of  the  predicted 
measurement  residuals  and  is  the  mechanism  for  adapting  the  filter  for 
maneuvering  targets.  The  transition  matrix,  F,  is  built  up  fnxu  sub- 
matrices, f: 


f 0 - 
0 f - 


where  the  submatrix  f is  given  by; 


I 0 T 0 


0 1 0 T 
0 0 1 0 


0 0 0 I 


and  T is  the  time  step;  t “ t[  - t^. 
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Measurements 


Four  kinds  of  measurements  can  be  used  in  MURLOC:  (1)  x-and-y 
position,  (2)  range -and-bearing,  (3)  bearing-only,  and  (4)  course-and- 
speed.  Standard  deviations  of  errors  in  the  measurements  are  also  input 
parameters.  Measurements  are  processed  each  time  step,  time  steps  may 
be  variable,  and  multiple  measurements  at  any  one  time  may  occur. 
Measurement  errors  are  assumed  to  be  uncorrelated  from  one  measurement 
to  another.  Also,  range  errors,  bearing  errors,  course  errors,  and  speed 
errors  are  assumed  to  be  uncorrelated  with  each  other.  However,  the 
x-and-y  position  errors  are  assumed  to  be  correlated  and  are  given  in 
terms  of  error  ellipse  parameters--the  standard  deviations  along  the  two 
principal  axes,  and  the  angle  from  North  to  the  major  axis.  Latitude- 
and-longitude  measurements  can  be  processed  by  first  transforming  to  a 
Cartesian  coordinate  system  a few  hundred  miles  in  size.  Latitude-and- 
longitude  measurements  are  the  assumed  origin  of  the  x-and-y  position 
measurements , 

Before  being  filtered,  the  measurements  are  transformed  to 
Cartesian  coordinates  and  the  coordinate  errors  are  approximated  by 
linear  functions  of  the  measurement  errors.  Usually  tracking  algorithms 
use  an  "extended"  Kalman  filter  that  linearizes  the  measurement  equations 
around  the  predicted  state  by  calculating  first-order  partial  derivatives 
This  procedure  is  acceptable  if  the  true  state  values  are  inside  the 
linear  region;  however,  if  the  predicted  state  is  in  error  by  a large 
amount,  the  tracker  will  behave  poorly.  Instead,  MURLOC  linearizes 
around  the  data  point  before  filtering,  and  thus  avoids  the  problem  of 
calculating  partial  derivatives  at  the  wrong  place. 

Range-and-bearing,  bearing-only,  and  course-and-speed  measure- 
ments are  transformed;  x-and-y  measurements  do  not  need  to  be  transformed 
For  example,  range-and-bearing  (r,0)  are  transformed  into  a measurement 
vector,  y,  by: 


y = 


r sin  9 
.r  cos  0 . 
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The  measurement  covariance,  R,  is  calculated  by  an  approximation  that 

assumes  that  the  standard  deviations  of  the  errors  in  range  and  bearing 

(o  ,a  ) are  small: 
r 


cos  0 sin  9 
-sin  0 cos  9 


r 2 2 „ 
r Og  0 


0 a 


r J 


cos  0 -sin  6 
sin  0 cos  6 


Bearing-only  measurements  require  an  estimate  of  range  before 
the  transformation  can  be  performed.  The  range  is  estimated  from  the 
predicted  positions  and  covariance  of  the  observer  and  target  units. 

The  estimated  range  is  calculated  by  expanding  (or  contracting)  the 
predicted  relative  error  ellipse  until  the  ellipse  just  touches  the 
measured  bearing  line;  the  tangent  point  then  defines  the  estimated 
range,  as  shown  in  Figure  2.  The  "predicted  relative  error  ellipse" 
is  the  error  ellipse  that  is  relative  to  the  observer;  it  is  computed 
from  the  predicted  covariance  elements  of  the  observer  and  target  (the 
equations  are  in  Subroutine  REPORT,  which  is  in  Appendix  A).  The 
standard  deviation  of  range  error  is  assumed  to  equal  the  estimated 
range.  Thus,  as  shown  in  Figure  2,  the  bearing-only  measurement  is 
transformed  to  a long,  thin  ellipse  that  lies  along  the  measured  bearing 


FIGURE  2 LINEARIZATION  OF  BEARING-ONLY  MEASUREMENT 
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and  is  centered  on  the  estimated  range.  The  measurement  vector,  y,  and 
the  covariance  matrix,  K,  are  then  computed  as  though  the  bearing-only 
measurement  were  a range-and-bear ing  measurement. 

Even  though  a psuedo-range  is  used  as  data,  and  the  assumption 
of  small  range  error  is  violated,  the  above  metliod  of  tracking  witli 
bearing-only  data  seemed  to  work  quite  well.  Another  method  for  estimating 
range  was  tried  but  it  introduced  large  range  errors  after  a target 
maneuver.  The  other  method  estimated  the  range  by  the  intersection  of  a 
line  that  was  perpendicular  to  the  bearing  line,  and  contained  the 
predicted  target  position.  The  other  method  was  computationally  simple, 
but  the  "ellipse  tangent"  algorithm  shown  in  Figure  2 proved  to  be  far 
superior  to  the  "perpendicular  intersection"  algorithm. 

Course-and-speed  (\,s)  measurements  are  transformed  in  the  same 
way  as  range-and-bear ing  measurements,  except  that  the  measurement  vector 
consists  of  velocity  components: 


y = 


s sin  Y 

s cos  Y 


The  covariance, 

R = 


R,  is  also  similar 

cos  Y sin  Y 
-sin  Y cos  Y 


to  the  range-and-bearing 


cos  \ -sin 
sin  > cos 


case : 

“ 

\ 


where  the  input  parameters  and  02  are  the  principal -axis  components 
of  the  velocity  error  elipse:  is  perpendicular  to  the  velocity 

vector,  and  O2  is  parallel  to  the  velocity  vector.  This  method  of  input 
is  used  instead  of  course  and  speed  sigmas  so  that  zero-velocity  error 
distributions  can  be  input  easily.  Course  and  speed  measurements  are 
thought  of  as  a polar-coordinate  representation  of  the  velocity  vector; 
the  error  along  the  vector  is  assumed  to  be  uncorrelated  with  the  error 


I*. 
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across  the  vector.  When  the  course  and  speed  standard  deviations 

(a  ,0  ) are  small,  the  approximation: 

> s 


s 0. 


*^2 


o 

s 


can  be  used  to  estimate  the  velocity  error  ellipse. 

By  linearising  the  input  data,  the  measurement  vector,  y, 
becanes  a random  vector  that  is  a linear  function  of  the  true  state 
values,  X: 


y = M X + V 

where  M is  the  measurement  matrix,  and  v is  a random  vector  that  is 
normally  distributed  with  zero  mean  and  covariance,  R, 

I’he  value  of  the  measurement  matrix,  M,  depends  on  the  kind 
of  measurement  and  the  units  involved,  M is  represented  as  a matrix  of 
submatrices: 


M = I ^ 


where  k “ l,2,,,,n  is  the  index  of  a unit  being  tracked,  f'or  x-and-y 
position  measurements,  M is  determined  by: 


A = 


10  0 0 
0 1 0 0. 


I A for  k = index  of  the  reporting  unit 
I 0 otherwise , 


For  range-and-bearlng  or  bearing-only  measurements,  M is  determined  by: 


M 


k 


I A for  k = index  of  the  target  unit 
\ -A  for  k = index  of  the  observing  unit 
\ 0 otherwise. 
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For  course-and-speed  measurements,  M is  determined  by: 


B = 


'0010' 

.0  0 0 1. 


I" 

|o 


for  k = index  of  the  reporting  unit 
otherwise . 


3 . Predicted  State 

MURLOC  assumes  that  the  units  move  with  constant  course  and 

speed.  The  estimated  state  at  time  t = tj^  is  predicted  from  the  estimated 

state  at  t : 

o 

X = F X 

o 

T 

The  predicted  covariance  is  not,  however,  the  usual  F P F matrix  with 

o 

some  constant  model  noise,  Q,  thrown  in.  Instead,  the  predicted  co- 
variance  is  required  to  be  consistent  with  the  measurement  by  making  Q 
a function  of  the  measurement  residuals.  The  residual  vector,  r,  is 
supposed  to  be  normally  distributed  with  zero  mean  and  covariance  matrix,  S 


r = y - M X 

S = R + M F P F^ 
o 

However,  when  a unit  changes  course  or  speed,  the  residual  vector  can  be 
many  standard  deviations  from  zero  and  the  statistics  of  r are  not 
what  they  are  supposed  to  be  because  the  assumptions  of  the  mod^i  are 

1* 

violated,  A procedure  that  is  philosophically  similar  to  Jazwinski's 
(but  not  mathematically  the  same)  is  used  to  let  the  model  noise  Q 
adapt  to  the  residuals. 

The  adaptive  algorithm  is  very  simple.  First  a beta  factor  is 
calculated: 

0=1-  exp  (-i  r^  S ^ r). 


References  are  listed  at  the  end  of  this  report. 
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An  error  vector,  q,  is  then  defined  in  terms  of  the  residual  vector, 
but  dimensioned  the  same  as  x: 

q = ( 0 0 ...  qj^  ...  00 

where  k is  the  index  of  the  target  unit.  The  vector  q^^  depends  on 
whether  the  measurement  is  a position  measurement  (x-and-y,  range-and- 
bearing,  or  bearing-only)  or  velocity  measurement; 

Position:  ^ ^ 

Velocity:  q^^  = ( 0 0 r^^  r^  ). 

Finally,  the  model  noise,  Q,  is  constructed  by  taking  the  outer  product 
of  q with  itself  and  reducing  the  resulting  matrix  by  the  beta  factor: 

Q = p q 

The  predicted  covariance  of  x is  then  given  by: 

P = F P F^  + Q . 
o 

The  purpose  of  the  adaptive  model-noise  matrix,  Q,  is  to  open 
up  the  predicted  covariance  so  that  the  most  recent  measurement  has  a 
significant  influence  on  the  filtered  state.  When  the  residual  vector 
is  large,  the  Q matrix  will  also  be  large,  P will  then  be  large,  and  the 
relatively  smaller  measurement  covariance,  R,  will  cause  the  filtered 
state  to  be  drawn  to  the  measurement.  This  idea  is  shown  in  Figure  3. 

Residuals  can  be  large  when  a unit  maneuvers;  unfortunately, 
residuals  can  also  be  large  for  large  but  natural  fluctuations  in  the 
measurement.  Thus,  there  is  a tradeoff  between  (1)  making  the  algorithm 
sensitive  to  target  maneuvers  and  letting  it  be  influenced  by  bad  data, 
and  (2)  decreasing  the  maneuver  sensitivity  so  that  bad  data  can  be 
smoothed  out.  For  the  beta  factor  as  defined  above,  the  algorithm  is 
rather  sensitive  to  the  residuals,  as  can  be  seen  in  Figure  4.  For  ex- 
ample, a two-sigma  residual  produces  B • 0.86  and  will  cause  a iiwderate- 
to- large  influence  on  the  predicted  covariance.  To  reduce  the  Influence 
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FIGURE  3 EXPANSION  OF  PREDICTED  COVARIANCE  IN  THE  DIRECTION  OF 
THE  RESIDUAL  VECTOR 
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FIGURE  4 MODEL  NOISE  COVARIANCE  REDUCTION  FACTOR 


of  residuals,  the  beta  factor  can  be  raised  to  some  power;  for  example, 

0 will  cause  the  algorithm  to  react  only  to  very  large  residuals. 

The  construction  of  the  Q-raatrix  was  Intuitive  rather  than 
based  on  a probabilistic  derivation.  Our  desire  was  to  create  an 
algorithm  that  was  continuous--that  is,  a Q-matrix  that  does  not  suddenly 
jump  in  value  as  the  residual  becomes  large. 

An  example  of  a Q-matrix  that  does  jump  is  the  following; 
let  Q equal  zero  if  the  residual  is  inside,  say, a two-sigma  elliptical 
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window;  and  let  Q equal  a given  noise  matrix  if  the  residual  is  outside 
the  window.  Another  similar,  but  more  sophisticated  method  is  Bayesian 
in  nature;  Let  the  algorithm  decide  between  a maneuvering  or  nonmaneuver- 
ing hypothesis.  The  maneuver  covariance  matrix,  the  a-priori  probabilities, 
and  a tactical  value-of-decision  matrix  can  be  input  to  reflect  the 
scenario  conditions.  In  both  of  the  above  methods,  the  algorithm  must 
choose  between  two  (or  more)  predicted  covariance  matrices  each  time  a 
measurement  is  taken.  These  decision-type  algorithms  can  underreact  to 
a slowly  maneuvering  target  that  stays  inside  the  window,  and  overreact 
to  bad  data  that  fall  outside  the  window.  In  an  attempt  to  smooth  out  the 
under-  and  overreaction  properties  of  the  decision-type  algorithm,  the 
continuous-Q  algorithm  was  created. 

T 

The  choice  of  q q as  the  basis  matrix  for  Q was  prompted 
by  the  desire  to  enlarge  the  predicted  state  covariance  only  in  the 
direction  that  is  indicated  by  the  data.  This  procedure  keeps  the  Q- 
matrix  as  small  as  possible  but  at  the  same  time  makes  the  predicted 
state  covariance,  P,  consistent  with  the  residual  vector,  r.  The 
velocity  components,  and  r^/T,  are  included  in  the  q vector 

because  the  residual  position  error  is  assumed  to  be  caused  by  a maneu- 
vering target  that  has  changed  course  and  speed;  therefore  the  velocity 
submatrix  of  the  predicted  covariance  must  be  increased  to  let  the 
state  velocity  move  to  a new  value. 

When  multiple  measurements  occur  at  the  same  time,  the  order 
in  which  they  are  filtered  affects  the  solution.  The  reason  is  that  the 
covariance  matrix,  P,  is  a function  of  the  data,  whereas  in  a nonadaptive 
Kalman  filter  the  covariance  is  independent  of  the  data,  and  the  order 
in  which  the  simultaneous  measurements  are  filtered  makes  no  difference. 
With  the  adaptive  algorithm,  the  solution  will  be  better  for  one  ordering 
of  simultaneous  measurements  versus  some  other  ordering.  The  best  order 
cannot  be  determined  as  the  calculations  are  made;  this  nonoptimum 
property  of  the  MURLOC  adaptive  algorithm  is  part  of  the  price  that  is 
paid  for  the  capability  of  tracking  a maneuvering  target  while  using  a 
constant -velocity  model. 


J 
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Two  ways  of  computing  simultaneous  measurements  were  tried. 

The  first  method  computed  the  Q-matrix  for  the  first  measurement  only, 
and  then  let  Q equal  zero  for  the  rest  of  the  measurements.  The  first 
method  did  not  perform  as  well  as  the  second  method,  in  which  Q was 
computed  for  each  of  the  measurements.  For  the  second  method,  the 
value  of  tau  in  the  velocity  error  computation  was  held  constant  and 
was  set  equal  to  the  time  duration  between  sets  of  measurements.  The 
second  method  caused  more  fluctuations  in  the  solution,  especially  in 
the  velocity  components,  but  it  seemed  to  react  to  a maneuvering  target 
more  quickly  than  the  first  method. 

4.  Filtered  State 

Once  the  predicted  state  vector  and  covariance,  x and  P,  were 
calculated,  standard  Kalman  filter  equations  were  used  to  compute  the 
filtered  state  vector  and  covariance  matrix,  Xj^  and  : 

K = P (R  + M P M^)'^ 

L=I  - KM  (1=  Identity  matrix) 

Xj^  = X + K r 

T T 

Pj^  = L P L + K R K . 

The  above  covariance  equation  was  used  instead  of  the  more  common 
Pj^  = (I  - K M)  P equation,  because  the  above  method  (1)  produced  a 
symmetric,  positive  definite  matrix,  and  (2)  was  less  sensitive  to 
numerical  errors  in  the  filter  gain  matrix,  K. 

B.  MURLOC  Example 

As  a way  of  describing  the  input  and  output  format  of  the  MURLOC 
computer  program,  a four-unit  scenario  was  created  as  shown  in  Figure  5. 
There  are  two  observer  ships,  one  helicopter,  and  one  target  ship.  The 
first  ship  (Unit  1)  measures  a sonar  bearing  to  the  target  ship  (Unit  4) 
and  launches  the  helicopter  (Unit  3)  to  relay  information  and  serve  as 
a radar  target.  The  first  ship  measures  target  bearings,  radar  range - 
and-bearing  to  the  helicopter,  and  own  course-and-speed.  The  second 


19 


100 


OWN  SHIP 
(Unit  1) 


20 


0 


0 


SISTER  SHIP 
(Unit  2) 


15  kt 


20  40  60  80  1 00 

X-POSITION  — nmi 


FIGURE  5 FOUR-UNIT  SCENARIO  FOR  THE  MURLOC  EXAMPLE 


ship  (Unit  2)  also  measures  target  bearings,  radar  range  and  bearing  to 
the  helicopter,  and  own  course  and  speed,  and  then  relays  this  informa- 
tion to  the  first  ship.  The  sonar  bearing  measurements  from  both  ships 
are  assumed  to  be  smoothed  for  10  minutes  and  have  a standard  deviation 
of  1.5  degrees.  Table  3 shows  the  assumed  sigmas  of  the  various 
measurements . 

Measurements  from  Unit  1 are  assumed  to  be  processed  5 minutes 
late  and  measurements  from  Unit  2 are  processed  10  minutes  late.  Table  4 
shows  the  time  sequence  of  measurements.  "Data  time"  is  the  time  when 
the  measurement  is  taken;  "filter  time"  is  the  time  when  the  measurement 
is  processed.  For  example,  at  time  t = 10,  Unit  2 measures  (1)  bearing 
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Table  3 


STANDARD  DEVIATION  OF  MEASUREMENT  ERRORS 
IN  THE  MURLOC  EXAMPLE 


Measurement 

Sigma 

Bearing  (1,4)* 

1.5  deg 

Bearing  (2,4) 

1.5  deg 

Range  (1,3) 

0.5  nmi 

Bearing  (1,3) 

2 deg 

Range  (2,3) 

0.5  nmi 

Bearing  (2,3) 

2 deg 

Course  (1) 

1.9  deg 

Speed  (1) 

1 kt 

Course  (2) 

3.8  deg 

Speed  (2) 

2 kt 

Bearing  (1,4)  is  the  bearing 
from  Unit  1 to  Unit  4,  etc. 


to  Unit  4,  (2)  range  and  bearing  to  Unit  3,  and  (3)  own  course  and  speed. 
These  data  are  relayed  to  Unit  1 via  the  helicopter,  and  10  minutes  later 
(t  = 20)  they  are  filtered  onboard  Unit  1.  Table  4 shows  that  the  example 
data  are  mostly  bearing  measurements  to  the  target  with  occasional  radar 
measurements  on  the  helicopter  and  course-and -speed  measurements  on  the 
ships. 

The  numbers  listed  in  Table  4 are  the  true  values  of  the  parameters 
calculated  at  the  data  time;  thus,  the  tracking  example  uses  measurements 
with  no  errors.  Tracking  results  produced  by  zero-error  measurements 
cannot  be  used  to  evaluate  the  tracking  algorithm,  but  instread  can  be 
used  to  roughly  estimate  the  magnitude  of  tracking  errors  that  would  be 
expected  if  the  measurements  were  randomized  around  the  true  values. 
Zero-error  measurements  were  used  for  the  MURLCXI  example  because  the 
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Table  4 


LIST  OF  MEASUREMENTS  FOR  THE  MURLOC  EXAMPLE 


T ime 
(min ) 

Sonar  Bearing 
(deg) 

Radar  Range  Bearing 
(nmi,  deg) 

Course 

(deg 

, Speed 
, kt  ) 

Data 

Filter 

Units  1,4 

Units  2,4 

Units  1,3 

Units  2,3 

Unit  1 

Unit  2 

0 

5 

081  .9 

000,  15 

5 

10 

083.6 

5.6,  122.2 

10 

20 

040.8 

54.1,  359.4 

000,  15 

15 

20 

087.2 

20 

30 

040.9 

25 

30 

091.4 

27.8,  122.2 

025,  18 

JO 

40 

041.2 

47.6,  016. 5| 

35 

40 

096.4 

1 

Effect  of  these  | 

40 

50 

041 .4 

measurements  shovm  in  . 

Appendix  B 

50 

100.1 

1 

50 

60 

040.0 

55 

60 

102.2 

60 

70 

038.5 

65 

70 

104.6 

70 

80 

036.9 

75 

80 

107.5 

80 

90 

038.5 

85 

90 

109.5 

90 

100 

043.3 

95 

100 

111.2 

100 

110 

048.0 

17.3,  126.1 

53.2,  014.7 

- - - . j 
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emphasis  of  this  section  is  on  the  input  and  outi)ut  fonnat  of  MlIRl.OC 
rather  than  its  performance.  Tracking  performance  is  the  subject  of 
the  next  chapter. 

1 . 

There  are  three  versions  of  the  MURl^OC  computer  program: 

(1)  a time-share  prototype  version  that  requires  the  data  to  be  entered 
interactively,  (2)  a generalized  version  that  uses  card  input,  and 
(3)  a specialized  version  that  uses  randomized  data  contained  on  a disc 
file.  The  card-input  version  of  MURLOC  is  listed  in  Appendix  A and  is 
the  subject  of  this  section.  The  file-output  version  of  MURLOC  was 
used  in  the  Monte  Carlo  performance  analysis  of  the  next  chapter. 

Table  5 shows  the  data  card  format  that  is  used  to  input 
measurements  to  the  program.  For  example,  the  card  for  the  bearing 
measurement  from  Unit  2 to  Unit  4 at  time  t = 10  reads: 


TP 

TD 

KEY 

KO 

K 

XI 

X2 

20.0 

10. 0 

3 

2 

4 

40.8 

1 .5 

Values  for  XJ  through  X5  do  not  need  to  be  punched  for  a bearing 
measurement  card. 

To  initialize  the  program.  Card  1 contained  TD  = 0,  KEY  = 6, 
and  KO  = 4.  This  card  told  the  program  to  initialize  the  state  at 
t = 0 and  to  set  the  number  of  units  to  4.  Qird  2 was  a position 
measurement  for  Unit  1 so  that  the  coordinate  system  could  be  defined. 
Card  3 was  a course -and-speed  measurement  on  Unit  1.  Card  4 Initialized 
the  position  of  the  target  by  using  a range-and-bear ing  fonnat;  bearing 
was  measured  (81.9'^)  and  range  was  estimated.  The  bearing-only  format 
could  not  be  used  for  the  first  bearing  measurement  because  tlie  position 
of  the  target  was  not  yet  defined;  therefore  a large,  but  reasonable, 
value  of  range  was  used-- in  this  case,  50  nmi,  as  compared  to  the  true 
value  of  70  nmi.  Card  5 set  the  course  and  speed  of  the  target  to  zero 
and  the  velocity  sigmas  to  30  knots.  Cards  6 through  34  contained  the 
data  shown  in  Table  4 starting  with  the  second  line.  The  final  card 
contained  KEY  ■ 7 to  end  the  run. 
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Table  5 

DATA  CARD  FXIRMAT 


Variable 

Meaning 

Units 

Coiumns 

TP 

Time  at  which  the  measurement 
is  filtered 

min 

1-10 

TD 

Time  of  measurement 

min 

11-20 

KEY 

Indicator: 

1 •=  Position  measurement 

2 = Range-and-bearing 

measurement 

3 = Bearing-only  measurement 

4 = Course-and-speed 

measurement 

5 = Output 

6 = New  run 

7 = End 

21 

KO 

Index  of  observing  unit 

- 

25 

K 

XI 

index  of  target 

Key: 

• 

27 

1 = x-position^ 

2 = Range 

3 = Bearing  (true) 

4 = Course 

nmi 

nnii 

deg 

deg 

31-40 

X2 

1 = y-position 

2 = Bearing  (true) 

3 = Bearing  sigma 

4 = Speed 

nmi 

deg 

deg 

kt 

41-50 

X3 

1 = Minor  axis  sigma 

2 = Range  sigma 

4 = Sigma  perpendicular  to 
velocity  vector 

nmi 

nmi 

51-60 

X4 

1 = Major  axis  sigma 

2 = Bearing  sigma 

4 “ Sigma  parailei  to 
velocity  vector 

nmi 

deg 

kt 

61-70 

X5 

1 = Angle  to  major  axis 
from  North 

deg 

71-80 

2.  Output 


A partial  listing  of  the  output  from  the  example  scenario  Is 
contained  In  Appendix  B.  The  program  prints  out  the  Information  con- 
tained on  a data  card  and  then  prints  the  results  of  filtering  those 
data.  The  filtered  results  are:  (1)  the  time  of  the  data  and  the  time 
of  the  calculation,  (2)  the  location  of  one  unit  relative  to  another, 
and  (3)  the  position  and  velocity  of  each  unit.  The  relative  location 
parameters  are: 

OBSR  UNIT  RNG  BRG  SIGl  SIG2  ANG  CEP. 

Under  the  OBSR  and  UNIT  headings  are  listed  all  combinations  of 
observer  and  target  unit  Indices.  RNG  and  BRG  are  the  "predicted  range 
and  bearing";  this  parameter  pair  Is  a polar  coordinate  representation 
of  the  vector  from  the  predicted  position  of  the  observer  to  the  target 
(RNG  and  BRG  are  not  average  values  of  range  and  bearing).  SIGl  and 
SIG2  are  the  predicted  minor  axis  and  major  axis  standard  deviations,  and 
ANG  Is  the  angle  from  North  to  the  major  axis  of  the  error  ellipse. 

The  error  ellipse  Is  In  a coordinate  system  that  Is  relative  to  the 
observer  unit.  CEP  Is  the  predicted  radius  of  a circle,  centered  on  the 
point  defined  by  the  predicted  range  and  bearing,  such  that  there  Is  a 
50  percent  probability  that  the  point  defined  by  true  values  of  range 
and  bearing  lies  Inside  the  circle.  CEP  Is  a single -parameter  measure 
of  the  three -parameter  error  ellipse,  and  Is  useful  In  quickly  assessing 
the  track  quality.  CEP  Is  a function  of  SIGl  and  SIG2;  the  equations 
are  given  In  Function  FCEP,  which  Is  In  Appendix  A. 

The  output  parameters  for  the  predicted  position  of  the  units 

are: 


UNIT  X Y SIGl  SIG2  ANG  CEP 


and  the  parameters  for  the  predicted  velocities  are: 

CRS  SPD  SIGl  SIG2  ANG  CEP. 
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The  predicted  x-and-y  position  of  each  unit  is  listed  along  with  the 
parameters  describing  the  position  error  ellipse  (SIGl,  SIG2,  ANG,  and 
CEP).  Error  ellipses  are  in  absolute  (geographic)  coordinates,  not  in 
relative  coordinates  as  are  the  range-and-bearing  position  error 
ellipses.  CRS  and  SPD  are  the  "predicted  course  and  speed"  of  each  unit; 
this  parameter  pair  is  a polar  coordinate  representation  of  the  predicted 
velocity  vector  (CRS  and  SPD  are  not  average  values  of  course  and  speed). 
The  next  four  parameters  (SIGl,  SIG2,  ANG,  and  CEP)  are  parameters  of  the 
velocity  error  ellipse.  A velocity  error  ellipse  is  used  because  this 
is  a more  natural  way  to  show  the  velocity  covariance  terms  than 
numerically  calculating  the  predicted  course  and  speed  error  standard 
deviations  and  correlation  between  course  and  speed  errors. 

Figure  6 shows  the  results  of  the  four-unit  scenario  example. 
The  range  from  the  first  ship  to  the  target  ship  is  plotted  against 
time.  The  target  maneuvers  at  t = 40  and  again  at  t = 75  minutes.  The 
first  maneuver  is  not  discernible  on  the  range  plot,  but  the  second 
maneuver  is  more  radical  and  is  easily  seen.  The  true  range  is  plotted 
and  the  tracking  results  using  zero-error  measurements  are  shown  as  a 
dotted  region.  The  region  shows  the  approximate  magnitude  of  position 
errors  that  would  be  expected  if  randomized  measurements  were  used.  The 
dotted  region  was  plotted  by  using  the  predicted  range  plus  and  minus 
the  relative-position  CEP  value  of  Unit  4 relative  to  Unit  1.  The 
stairstep  effect  is  caused  by  the  ever-expanding  predicted-position  CEP 
being  periodically  contracted  by  measurements  that  are  filtered  at 
10-minute  intervals.  The  range  estimate  overshoots  the  second  maneuver 
but  comes  back  and  brackets  the  true  range.  The  overshoot  is  mostly  due 
to  the  10-minute  delay  in  information  from  the  second  ship,  rather  than 
slow  reaction  of  the  tracker.  Figure  6 shows  that  position  errors  of 
roughly  5 nmi  can  be  expected  in  a four-unit  targeting  scenario  involving 
over-the-horizon  ranges  and  measurement  errors  as  shown  in  Table  2. 


TARGET  RANGE  — nautical  miles 


Ill  MURLOC  TRACKING  CAPABILITY 


MURLOC s tracking  capability  was  investigated  by  a Monte  Carlo 
analysis  of  tracking  errors,  A three-unit  scenario  was  defined  and 
50  replications  of  randomized  measurements  were  calculated  and  stored 
on  a disc  file.  Then  MURLOC  was  run  50  times  using  the  random  measure- 
ments; the  resulting  position  and  velocity  errors  were  averaged  and 
output  as  measures  of  tracking  capability.  In  addition  to  running  the 
MURLOC  algorithm,  two  versions  of  another  tracking  algorithm  were  run 
to  compare  against  MURLOC.  Thus,  by  using  the  same  file  of  50  replica- 
tions in  all  three  algorithms,  a valid  comparison  could  be  made.  The 
comparison  algorithms  employed  a batch-processed  least-squares  methodology 
that  is  described  in  part  in  this  chapter  and  in  detail  in  the  next 
chapter. 


A.  Three-Unit  Scenario 

Figure  7 shows  the  three-unit  scenario  used  for  the  tracking 
capability  analysis.  Two  ships  (Units  1 and  2)  cross-fix  the  target 
ship  (Unit  3)  with  bearing-only  measurements,  and  find  their  baseline 
with  range-and-bearing  measurements.  Measurement  sets  were  10  minutes 
apart;  there  were  10  sets  of  measurements  (10  time  points)  each 
replication,  and  8 measurements  were  made  within  each  measurement  set. 
Table  6 lists  the  8 measurements  that  were  generated  at  each  time  point, 
along  with  the  standard  deviations  used  in  the  randomizing  the  measure- 
ments, A normal  distribution  of  measurement  error  was  used  and  the  mean 
value  of  the  measurement  of  a parameter  was  set  equal  to  the  true  value 
of  the  parameter  at  the  time  point.  There  was  no  delay  time  between 
measurement  and  calculation;  in  the  vocabulary  of  the  previous  chapter, 
"data  time”  equaled  "filter  time."  The  two  observing  ships  maintained 
a constant  course  and  speed  throughout  the  scenario.  The  target  ship 
made  a 90-degree  course  change  at  the  fifth  time  step  (t  = 40  min),  but 
maintained  its  speed  of  24  knots. 
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FIRST  SHIP 
(Unit  II 


TIME  STEP  ■=  10  min 


FIGURE  7 THREE-UNIT  SCENARIO  FOR  THE  MURLOC  TRACKING  CAPABILITY  ANALYSIS 

At  each  time  step,  MURLOC  processed  the  8 measurements  as  5 groups 
of  simultaneous  measurements.  The  order  in  which  the  measurement  groups 
were  filtered  is  shown  in  Table  7.  Switching  the  order  of  Bearing 
(2,3)  and  Bearing  (1,3)  degraded  the  tracking  result  slightly.  Other 
than  investigating  the  above  bearing  switch,  no  attempt  was  made  to 
find  the  best  order  for  the  simultaneous  measurements. 

t 

I 

B.  Least-Squares  Comparison 

A least-squares  methodology  was  used  to  compare  against  MURLOC  by 
1 formulating  two  algorithms;  An  algorithm  that  decided  the  best  set  of 

I measurements  using  an  optimization  procedure,  and  (2)  an  algorithm  that 

[ was  told  when  the  target  maneuvered  and  could  therefore  process  the 

I best  set  of  measurements.  These  two  algorithms  are  called  "self- 

I optimizing"  and  "maneuver-known,"  and  are  discussed  below. 
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Table  6 


i 

i 

1 


STANDARD  DEVIATION  OF  MEASUREMENTS 
IN  THE  THREE-UNIT  SCENARIO 


Measurement 

Sigma 

1. 

Bearing  (1,3)* 

1.5 

deg 

2. 

Bearing  (2,3) 

1.5 

deg 

3. 

Range  (1,2) 

0.5 

nmi 

4. 

Bearing  (1,2) 

2 

deg 

5. 

Course  (1) 

I 

deg 

6 . 

Speed  ( 1 ) 

1 

knot 

7. 

Course  (2) 

1 

deg 

8. 

Speed  (2) 

1 

knot 

Bearing  (1,3)  is  the  bearing 
from  Unit  I to  Unit  3,  etc. 


Table  7 

ORDER  OF  FILTERING  THE  MEASUREMENTS  IN  MURLOC 


Order 


1 

2 

3 

4 

5 


Measurement  Group 


Bearing  (2,3) 

Range -and-bearing  (1,2) 
Bearing  (1,3) 

Course -and-speed  (1) 
Course-and-speed  (2) 


I 

i 

I 

i 

I 


I 
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1 . Self-Optimizing  Algorithm 

The  self-optimizing  least-squares  algorithm  processed  a 
variable  number  of  measurements  at  each  time  step.  For  example,  at  the 
n-th  time  step,  the  algorithm  first  processed  16  measurements: 

8 measurements  from  the  current  time,  t = n,  and  8 from  the  irrmediately 
preceding  time,  t = n - 1.  These  16  measurements  were  used  to  estimate 
10  state  parameters:  2 velocity  components  for  Unit  1,  4 position  and 
velocity  components  for  Unit  2,  and  4 components  for  Unit  3.  With  a 
state  estimate  based  on  16  measurements,  the  algorithm  then  calculated 
a test  statistic  and  saved  it.  The  process  was  repeated  using  24 
measurements  from  time  steps  t = n,  n - 1,  and  n - 2,  and  again  a test 
statistic  was  calculated  and  saved.  The  calculations  were  repeated  as 
many  times  as  it  took  to  use  up  all  of  the  measurements;  thus  at  the  n-th 
time  step,  n - 1 state  estimates  and  test  statistics  were  calculated. 
Finally,  the  measurement  set  that  produced  the  minimum  test  statistic 
was  picked  as  the  "best"  set  to  use  at  time  t = n. 

The  above  optimization  procedure  was  used  because  of  the 
maneuvering-target  problem.  If  tracks  were  always  straight  lines,  then 
the  '■alculation  of  state  estimates  would  utilize  as  many  measurements 
as  possible.  But  the  chance  of  a maneuvering  target  required  the 
algorithm  to  decide  what  measurements  to  use,  since  measurements  taken 
before  a maneuver  cause  large  errors  in  estimating  the  state  after  the 
maneuver.  Ideally,  the  optimization  procedure  would  process  all 
measurements  from  the  time  of  the  maneuver  onward.  However,  in  practice, 
the  test  that  is  used  to  determine  the  measurement  set  necessarily 
operates  on  random  values  and  cannot  always  determine  the  truly  best  set 
of  measurements. 

2.  Maneuver -Known  Algorithm 

The  maneuver -known,  least-squares  algorithm  was  allowed  to 
process  the  truly  best  set  of  measurements  each  time  step.  For  example, 
at  Time  Step  9,  the  measurements  at  9,  8,  7,  6,  and  5 were  used,  but 
measurements  before  Step  5 were  not  used  because  they  occurred  before  the 
target  maneuver.  The  reason  for  defining  the  maneuver -known  algorithm 
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was  to  show  the  best  possible  use  of  the  measurements.  Since  it  was  a 
batch-processed  least-squares  methodology  that  is  allowed  to  know  the 
time  of  target  maneuver,  the  algorithm  produced  a lower  bound  on  the 
tracking  error.  It  must  be  remembered,  though,  that  MURLOC  and  the 
self-optimizing  least-squares  algorithm  were  real  trackers,  whereas  the 
maneuver-known  algorithm  was  not  real,  in  the  sense  that  it  could  not 
calculate  state  estimates  based  solely  on  measurements. 

C.  Results 

The  position  and  velocity  errors  as  a function  of  time  were 
averaged  over  50  replications  and  used  to  determine  tracking  capability. 

1 . Position  Errors 

"Position  error"  is  defined  as  the  distance  from  the  true 
position  of  the  target  to  the  estimated  position  in  a coordinate  system 
that  is  relative  to  Unit  1.  Position  error  is  thus  relative  position 
error,  and  it  is  the  magnitude  of  the  vector  difference  of  the  estimated 
position  vector  minus  the  true  position  vector,  using  Unit  1 as  the 
origin  for  the  two  vectors. 

Figure  8 shows  the  average  position  error  as  a function  of  time 
for  MURLOC  and  the  two  least-squares  algorithms.  MURLOC  performed  very 
well  when  compared  to  the  lower  bound. 

The  fact  that  MURLOC  and  the  self-optimizing  algorithm 
performed  about  the  same  suggests  that  the  optimizing  procedure  could  be 
improved  for  the  least-squares  tracker.  The  test  statistic  that  was 
used  is  related  to  the  sum-of-squares  of  the  residuals.  Another  test 
statistic  is  suggested  in  the  next  chapter,  but  we  did  not  have  the 
resources  to  investigate  it  or  other  alternative  self-optimizing 
procedures . 

Figure  9 shows  the  percent  of  replications  for  which  the  posi- 
tion error  is  less  than  the  estimated  CEP  of  Unit  3 relative  to  Unit  1. 
For  each  replication  and  at  each  time  step,  CEP  is  estimated  from  the 
covariance  matrix.  This  single-parameter  estimate  of  position  error  is 
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FIGURE  8 AVERAGE  POSITION  ERROR  COMPARISON  OF  MURLOC  TO  LEAST-SQUARES 
ALGORITHMS 


useful  in  making  targeting  decisions;  however,  CEP  is  a random  variable 
and  may  not  represent  the  actual  position  errors.  Figure  9 indicates 
that,  in  both  MURLOC  and  the  self-optimizing  least-squares  algorithm, 
the  CEP  values  were  smaller  than  they  should  have  been.  If  CEP  were 
estimated  better,  approximately  50  percent  of  the  replications  would 
result  in  position  errors  less  than  the  CEP  values. 


The  reason  for  the  overly  optimistic  estimate  of  the  error  is 
not  known,  but  since  the  maneuver-known  algorithm  approached  the 
theoretical  curve  better  than  the  other  two  algorithms,  optimun  use  of 
the  data  may  play  an  important  part  in  correctly  estimating  error 
variances . 


FIGURE  9 POSITION  CEP  PREDICTIVE  CAPABILITY 


2,  Velocity  Errors 

"Velocity  error"  is  the  distance,  in  velocity  space,  from  the 
true  velocity  point  to  the  estimated  velocity  point  of  Unit  3.  In  other 
words,  velocity  error  is  the  magnitude  of  the  vector  difference  of  the 
estimated  velocity  vector  and  the  true  velocity  vector  of  the  target. 

Figure  10  shows  the  average  velocity  errors  as  a function  of 
time  for  the  3 tracking  algorithms.  Even  at  best,  10  knots  error  is  not 
too  impressive,  for  MURLOC  or  the  self -optimizing  algorithm,  but  the 
long-range  geometry  combined  with  the  1.5-degree  bearing  sigmas  make  the 
estimation  of  velocity  most  difficult.  The  maneuver -known  algorithm 
results  are  significantly  lower,  and  show  the  improvement  that  might  be 
possible  with  an  accurate  maneuver  detector. 

Even  though  the  maneuver  is  known  to  occur  at  t = 40  min,  a 
large  velocity  error  is  recorded  for  the  maneuver -known  algorithm 
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FtGURE  W AVERAGE  VELOCITY  ERROR  COMPARISON  OF  MURLOC  TO  LEAST-SQUARES 
ALGORITHMS 


because  the  velocity  estimate  was  based  on  data  from  Time  Steps  1 
through  5,  The  target  course  estimate  was  close  to  315  degrees,  but  the 
program  defined  the  target's  course  at  Step  5 as  045  degrees,  whereas 
actually  the  course  is  both  315  and  045  degrees.  Thus,  a large  velocity 
error  shown  at  the  time  of  maneuver  is  a result  of  a programming  defini- 
tion, rather  than  a result  of  estimation  difficulties.  These  comments 
also  apply  to  the  two  other  algorithms. 

Figure  11  shows  the  effect  on  the  MURLOC  results  when  the 
bearing  errors  for  both  Bearing  (1,3)  and  Bearing  (2,3)  are  drawn  from 
a normal  distribution  with  a standard  deviation  of  1 degree  instead  of 
1.5  degrees.  The  improvement  indicates  that  the  velocity  errors  are 
significantly  affected  by  the  underlying  error  structure  of  the  scenario. 

Figure  12  shows  the  percent  of  replications  that  had  velocity 
errors  less  than  the  estimated  velocity  CEP.  The  results  show  that 
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FIGURE  11  EFFECT  OF  SIGMA  BEARING  ON  MURLOC  VELOCITY  ERROR 


MURLOC  and  the  self-optimizing  algorithm  were  as  overly  optimistic  in 
velocity  estimates  as  they  were  in  position  estimates.  The  maneuver- 
known  algorithm  again  approximates  the  theoretical. 

3 . Summary 

The  Monte  Carlo  results  for  the  defined  scenario  show  that 
MURLOC  can  track  a maneuvering  target,  that  the  target  position  errors 
in  relative  coordinates  will  be  close  to  the  lower  bound,  and  that  the 
target  velocity  errors  will  be  reasonably  small,  as  compared  to  a 
self-optimizing  least-squares  method.  MURLOC  estimates  its  own  errors 
smaller  than  they  really  are,  and  therefore  decisions  based  on  error 
estimates  would  have  to  take  MURLOC s overly  optimistic  behavior  into 
account.  An  examination  of  other  scenarios  and  other  least-squares 
optimization  procedures  would  undoubtedly  add  to  the  knowledge  of 
MURLOC s tracking  capability,  but  the  above  basic  analysis  does  show 
that  MURLOC  is  a reasonably  good  tracking  algorithm. 

37 


TIONS  WITH  VELOCITY  E 


IV  THE  COMPARISON  ALGORITHM 


The  tracking  capability  of  MURLOC  was  compared  to  the  capability  of 
another  tracking  algorithm  that  employed  a batch-processed  self-optimizing, 
least-squares  methodology.  Upon  receiving  measurements  at  the  n-th  time 
step,  the  self-optimizing  algorithm  computed  an  estimate  of  the  state 
vector  using  measurements  from  time  n and  n - 1;  then  it  computed  another 
state  estimate  using  measurements  from  time  n,  n - 1,  and  n - 2;  and  so 
on  down  to  the  first  time  step.  The  n - 1 state  estimates  were  checked 
to  see  which  one  produced  the  smallest  value  of  a defined  test  statistic; 
the  chosen  one  was  then  used  as  the  best  state  estimate  at  time  step  n. 
State  estimates  were  calculated  by  processing  measurements  in  a nonlinear 
weighted  least-squares  algorithm.  The  algorithm  employed  an  iterative 
method  to  find  the  best  least-squares  solution.  The  sections  that  follow 
describe  the  comparison  algorithm  in  detail. 


A.  Least-Squares  Problem 

At  each  time  step,  n,  a state  estimate,  z 


pn 


(where  p=lton-l). 


is  computed  such  that  a sum  of  weighted  squares  is  minimized.  The  sum 
of  squares,  is  given  by: 


G 

pn 


for  time-step  summation  indices  k ^ p to  n,  and  measurement  indices 
1 ■ 1 to  8.  The  actual  state  vector  at  time  n is  denoted  z , the 

j " 

measurements  at  time  k are  denoted  mj^,  and  the  measurement  models  at 
time  k are  denoted  ; these  three  objects  are  defined  below.  The 

measurement  variances  at  time  k are  denoted  Mj";  they  are  input  parameters 


There  are  n - 1 state  estimates,  z , at  time  step 

' pn'  ^ 


n,  and  one  of 


them  is  better  than  the  others.  The  best  one  is  found  by  computing  a test 


39 


staClsCic, 

value  S . 
pn 


S 


, and  choosing  the  state  estimate  that  produces  the  smallest 
The  test  statistic  is  defined  in  Section  C of  this  chapter. 


1.  State  Vector 

At  each  time,  n,  where  n ranges  from  2 to  N,  an  estimate  of 
position  and  velocity  for  each  of  3 units  is  required.  This  estimate  is 
based  on  a contiguous  time-sequence  of  measurements.  Let 

X . = The  x-coordinate  of  Unit  i at  time  n 
ni 

y . = The  y-coordinate  of  Unit  i at  time  n 
ni 

• 

X . = The  velocity  in  the  x-direction  at  time  n 
ni 

y^^  = The  velocity  in  the  y-direction  at  time  n 

where  i = 1,2,3  is  the  unit  index.  Position  is  expressed  in  nautical 
miles  and  velocity  is  expressed  in  knots.  The  state  vector  has  10  com- 
ponents and  is  defined  as: 


= (X 


nl  ^nl  ^2  yn2  *n2  yn2  ^n3  yn3  *n3  yn3>^ 


It  is  convenient  to  let  z . denote  the  J-th  component  of  z . 

nj  n 

The  state  vector  does  not  contain  the  x and  y positions  of 
Unit  1 because  only  the  relative  positions  need  to  be  estimated.  In 
effect,  the  x, y-coordinate  system  is  redefined  at  each  time  step  rela- 
tive to  the  position  of  Unit  1. 


2.  Measurement  Vector 

The  measurement  vector  at  time  n is  defined  as: 


m 

n 


n/ 


where  the  m^  (j  * 1 to  8)  are  defined  in  Table  8.  We  assume  that  m is 
n n 

a random  vector  distributed  as  a multivariate  normal  deviate  and  that 
the  components  of  m are  uncorrelated.  We  call  M the  covariance  matrix 

i " 

of  m and  use  M to  denote  the  j-th  diagonal  term  of  M . 
n n j a 
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Table  8 


DEFINITION  OF  MEASUREMENTS  FOR  THE  COMPARISON  ALGORITHM 


Symbol 

Meaning 

1 

★ 

"V 

Bearing  (1,2) 

2 

Bearing  (1,3) 

3 

Course  (1) 

4 

"he 

Range  (1,2) 

5 

"he 

Speed  (1) 

6 

"he 

Bearing  (2,3) 

7 

"he 

Course  (2) 

8 

"he 

Speed  (2) 

*Bearing  (1,2)  is  the  bearing  from 
Unit  1 to  Unit  2 at  time  step  k. 


3.  Measurement  Model 


A model  of  the  measurements  m^  at  time  k can  be  constructed 
from  the  state  vector  and  the  position  of  Unit  1 at  time  n . First 
define  the  relative  x-position  by: 


- X ) + (k  - n)  T (X  - X ) 
nj  nl  nj  ni 


where  t is  the  time-step  duration,  and  i and  j are  unit  indices.  Define 

nlc 

a similar  equation  for  y . Then  components  of  the  measurement  model  can 

1 8 

be  written;  for  example,  the  components  corresponding  to  m^^  and  are: 


F^(Zn)  = tan'^ 

+ (yn2>^f  • 


We  say  that 


J 

F,  (z  ) = (f!-(z  ) ...  F®(z  )) 
k n \ k n k n / 


is  the  measurement  model  at  time  k. 

B.  Nonlinear  Least-Squares  Algorithm 

References  that  surveyed  literature  directly  addressing  the 
solution  of  nonlinear- least  squares  problems  were  examined;  see  Bard,^ 
Broyden,^  Draper  and  Smith, ^ Fletcher,^  and  Powell.^  The  examination 
revealed  that  all  algorithms  that  solve  nonlinear  least-squares  problems 
are  iterative.  They  calculate  a sequence  of  points  Vj^,  V2,  ...  , that 
should  converge  to  a V*,  which  solves  the  nonlinear  least-squares 
problem.  These  algorithms  frequently  substitute  linear  approximations 
evaluated  at  for  the  nonlinear  functions,  and  these  linear  approxima- 
tions are  then  used  to  calculate  V,  , . A conclusion  of  this  review  of 

k+1 

relevant  literature  was  that  the  Gauss-Newton  method,  and  variants  of 

it,  were  the  best  algorithms  with  which  to  minimize  the  sum  of  squares, 

G . The  particular  variant  of  the  Gauss-Newton  method  that  we  used  in 
pn 

the  optimized.  Iterative,  least-squares  tracking  algorithm  makes  use  of 

a technique  to  ensure  that  G strictly  decreases  for  each  point  of  the 

pn 

sequence  generated. 

In  the  following  sections,  G is  denoted  G(z) ; in  other  words, 

pn 

the  p and  n subscripts  are  suppressed  and  the  functional  dependence 
on  the  state  vector  z s z^  is  emphasized. 
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1 


Gauss-Newton  Method 


1. 


ALgorithm  (a)  in  Bard^  was  used  as  the  variant  of  the  Gauss- 
Newton  method  for  handling  the  nonlinear  problem.  For  the  first  itera- 
tion, the  algorithm  was  provided  an  initial  estimate  of  the  state.  To 
reduce  the  chance  of  a nonconverging  solution,  we  used  the  true  state 
values  to  start  the  Gauss-Newton  algorithm. 

At  time  step  n and  iteration  i , is  approximated  by 

the  linear  function; 


fj(z)  = Fj[z(i)]  +'£  J^g(i)  [z  - z(i)] 

where  the  state  summation  index  is  s = 1 to  10,  the  measurement  index  is 
j * 1 to  8,  the  time  step  index  is  k = p to  n,  and; 


■ -aH- 


ns 


z * z (i) 
n n 


The  values  used  to  approximate  G(z)  by  the  sum  of  squares: 


k j 

The  vector  d(i),  which  minimizes  g(z),  is  then  determined. 

The  estimate  of  state  at  iteration  i is  denoted  z(i).  A new 
estimate  z(i+l)  is  found  using  the  formula: 


z(i+n  = z(i)  + r d(i) 


where  d(l)  is  the  solution  of  least-squares  problem  at  time  n and 
iteration  i,  and  r is  a scalar  chosen  so  that; 

G[z(i+1)]  < G[  z(i)]  . 


I 

j 

i 

I 


i 


1 

i 

j 

ji' 


ii- 

i 
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The  iterations  at  time  n continue  until  any  one  of  the  following  four 
termination  criteria  is  satisfied; 

• The  number  of  interpolations  (defined  below)  reaches  its 
maximum  allowed  value  (10  were  allowed)  and 

(J[z(i)  + r d(i)]  G[z(i)] 

• The  iteration  counter,  i,  readies  its  maximum  allowed 
value  (5  iterations  were  allowed). 

• All  absolute  components  |r  d(i)|  are  less  than  0.0001 
times  the  corresponding  absolute  components  |z(i)|  + 0.001; 
see  Bard,^  page  168. 

• G[z(i)]  - 0.001  < Gl  z(i+l)l  <G[z(i)] 

The  scalar  r at  iteration  i is  calculated  by  the  following  steps: 

(1)  Use  the  value  of  the  interpolation  index  q from 
iteration  i - I (if  i = 1 set  q = 0). 

(2)  Divide  q by  2 and  retain  the  integral  part. 

(3)  Set  r = 2"*'. 

(4)  If  G[z(i)  + r d(l)]  <Glz(i)],  use  z(i+l)  = z(i)  + 
r d(i)  as  the  new  estimate  of  state.  Otherwise 
perform  Steps  5 through  8. 

(5)  Construct  as  a function  of  s the  polynomial  of 
degree  2 that  equals  Gl  z(i)  + s d(i)]  for  the  two 
points  s ||  0 and  s ■ r,  and  for  which  the  slope 
equals  h V h at  s 0.  The  vector  li  and 

tlie  matrix  V are  defined  in  the  next  section. 

(6)  Compute  tlie  value  s at  whicli  the  first  deriva- 
tive of  the  polynomial  is  zero  and  call  this 
value  Sj^.  Let 

s^  = max  [0. 25r,  min  (0.75  r,  s^] 

An  execution  of  Step  6 is  called  an  "interpolation." 

(7)  Increase  q by  one. 

(8)  Replace  r with  S2  and  return  to  Step  4. 

2.  Least-Squares  Solution 

Several  algorithms  were  available  to  compute  d(i),  (see 
Lawson  and  Hanson).^  Since  the  covariance  matrix  of  d(i)  was  desired, 
an  algorithm  was  selected  that  provides  the  covariance  matrix  as  a 
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byproduct  of  the  computation  of  d(i).  Execution  time  and  programming 
ease  were  also  considered.  Consequently,  the  normal  equations  of  the 
linear  weighted  least-squares  problem  were  derived,  the  Cholesky  (or 
square-root)  method  was  applied  to  the  matrix  determined  by  the  normal 
equations,  the  upper  triangular  Cholesky  factor  was  inverted,  and  the 
d(i)  that  solved  the  linear  weighted  least-squares  problem  and  its 
covariance  matrix  was  then  computed. 


Below  we  show  the  form  of  the  normal  equations  for  the  linear  least- 
squares  problem.  Let  p and  n be  integers  such  that  p = 1 to  n - 1. 
Let  the  8-by-lO  matrix,  J (i),  contain  the  (j,s)-th  component, 

j ^ 

Jj^g(i)  for  k = p to  n.  Define  the  matrixes: 


Define  the  vector: 


- [/(«  ... 


T 

Jp(i) 


M 


M 


’M 


mn  - Fj2(i)] 


m - F fz(i)] 
P P 


The  normal  equations  at  time  n and  Iteration  i are  then  given  by: 

T -1  T -1 

J M J d(i)  = M D . 

T -1 

Now  J M J is  symmetric  and  positive  definite;  thus  Cholesky* s 
method  (see  Lawson  and  Hanson^  or  Forsythe  and  Moler®)  is  used  to  find 
an  upper  triangular  matrix  U such  that: 

T T -1 

U U - M J 
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aad  U is  invertible  by  using  formulas  from  Lawson  and  Hanson.  The 
covariance  matrix  of  the  solution  of  the  linear  weighted  least-squares 
problem,  V,  is  computed  from  the  formula: 

-1  -IT 
V = U (U  ) 

The  vector  h is: 

T -1 

h = J M D 

and  the  solution  of  the  linear  weighted  least-squares  problem,  d(i), 
is  calculated  by: 


d(i)  = V h . 


C.  Test  Statistic 

Two  test  statistics  are  presented  below;  the  first  was  used  in  the 
self-optimizing  least-squares  algorithm  presented  in  Chapter  III.  The 
second  statistic  is  suggested  as  a possible  candidate  for  the  optimiza- 
tion process,  but  no  experimentation  was  performed  to  see  if  it  would 
work  better  than  the  first  statistic. 


1.  Residual  Mean  Square 

The  residual  mean  square,  S,  was  the  test  statistic  minimized 

by  the  self-optimizing  algorithm.  For  each  time  step  n,  the  statistic 

Spn  (where  p = 1 to  n - 1)  was  computed  and  the  estimate  of  state, 

which  produced  the  smallest  S was  chosen  as  the  best  estimate  of 

pn 

state  at  time  n.  The  statistic  was  calculated  by: 

V “ - p + 1)  - 10] 


where  the  sum  of  squares, 

estimates,  z 

’ pn 


was  evaluated  at  each  of  the  state 


2.  Multiple  Correlation  Coefficient 

2 

The  square  of  the  multiple  correlation  coefficient,  R , Is  a 

statistic  that  Is  frequently  computed  for  multiple  regression  analysis 

2 

(see  Draper  and  Smith  ).  R Is  a measure  of  how  well  the  regression 

2 

explains  the  data:  R Is  between  0 and  1,  and  the  larger  It  Is,  the 

2 

better  the  regression.  R Is  offered  here  as  a candidate  test  statistic 

to  determine  the  best  of  the  n - 1 state  estimates,  z , Suppose  z 

’ pn  pn 

was  obtained  on  Iteration  1;  then  the  test  statistic  Is  given  by: 

2 T / T -1 

R - (h^  V h - a)  / (D  M D - a) 


where 


8(n  - p + 1) 


for  time-step  Indices  k - p to  n and  measurement  Indices  j = 1 to  8. 

2 

The  state  estimate  that  Is  associated  with  the  largest  value  of  R 

pn 

(where  p “ 1 to  n - 1)  would  then  be  chosen  as  the  best  estimate  at 
time  n. 
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Appendix  A 

MURLOC  COMPUTER  PROGRAM 


PROGRAM  MURLOC ( 1 NPUT , OUTPUT , TAPES* I NPUT , TAPE6=OUTPUT ) 

C 

C ««  KEY 

C ««  1 GEOGRAPHICAL  POSITION  DATA 
C ««  2 RANGE  AND  BEARING  DATA 

C ««  3 BEARING-ONLY  DATA 
C ««  4 COURSE  AND  SPEED  DATA 

C ««  5 OUTPUT 

C ««  6 NEW  RUN 

C *«  7 END 

C ««  K UNIT  NUMBER 

C ««  KO  OBSERVER  NUMBER 

C ««  TM  UPDATE  TIME  (MIN) 

C ««  RNG  RANGE  (NMI) 

C *«  BRG  TRUE  BEARING  (DEG) 

C ««  RSIG  RANGE  STANDARD  DEVIATION  (NMI) 

C «<  BSIG  BEARING  STANDARD  DEVIATION  (NMI) 

C ««  XE  EAST  COORDINATE  (NMI) 

C ««  YE  NORTH  COORDINATE  (NMI) 

C *«  R1  X- ROTATED  PRINCIPLE  STANDARD  DEVIATION  (NMI) 

C ««  R2  Y -ROTATED  PRINCIPLE  STANDARD  DEVIATION  (NMI) 

C «>  ANG  ANGLE  OP  ROTATION  FROM  NORTH  (DEG) 

C ««  CRS  COURSE  (DEG) 

C ««  SPD  SPEED  (KT) 

C «»  TO  LAST  TIME  OF  UPDATE  (MIN) 

C »«  Ed)  STATE  VECTOR  ( XI , Y1 , VX1 , VY1 . X2,  Y2.  VX2,  VY2,  ETC.) 

C *«  V(dJ)  COVARIANCE  MATRIX 

C »«  H(I,J)  STATE-TO-DATA  TRANSFORMATION  MATRIX 

C ««  D(N)  DIFFERENCE  OF  DATA  AND  PREDICTED  DATA 

C »«  W(I,J)  DATA  COVARIANCE  MATRIX 

C 
C 

COMMON/A/  KEY, TM, TO,KO.K,KM, IMAX,UB,TDO 
COMMON/B/  D(2),W(2,2),H(2,36) 

COMMON/ C/  E(36),EE(36),V(36,36),VV(36,36),A(36,36) 

COMMON/E/  RNG , BRG , RS I G, BS I G, XE , YE 
DIMENSION  X(5) 

C 


oo  oo  oo  oo  oo 


I 


««  GEOGRAPHICAL  DATA 
30  CONTINUE 

IF(TM.GT.TO)  CALL  UPDATE! 1 ) 

XEsXd)  SYE^XCZ) 

R1sX(3>  SR2=X(4)  SAN0cX(9) 

123K«4-2  si  1 >12*1 
0( 1 )sXE-E( I 1 ) 

0(2)>YE'E( 12) 

IM2=2*IMAX 
00  33  I < 1 , I M2 
33  H( I )s0. 

H(1,M)  = 1.  SH(2,I2)  = 1. 

B=ANGxUB 

CALL  ROTATE! R1,R2,B.W( 1,1 ) , W! 1 , 2) , W! 2, 2) ) 
W!2, 1 )xW!1 ,2) 

CALL  KALMAN! 2, F) 

LL*2  STM=TP  SGO  TO  110 

**  RANGE  AND  BEARING  DATA 
40  CONTINUE 

RN0xX!1)  SRSI03X!3) 

BR0«X!2)  SBSI0xX!4) 

CALL  REPORT 
CALL  KALMAN ! 2, F) 

LL«2  STMxTP  SGO  TO  1 1 0 

««  BEARING  ONLY  DATA 
44  CONTINUE 
BRGxX! 1 ) 

BSI0>X!2) 

CALL  REPORT 
CALL  KALMAN! 2, F) 

LL»2  STMxTP  SGO  TO  110 

»*  VELOCITY  DATA 
50  CONTINUE 

CRS-X!1)  SSP0xX!2) 

S1«X!3)  SS2>X!4) 

IFITM.GT.TO)  CALL  UPDATE! 1 ) 

B>CRS»UB 

I4xK*4  SI3xI4-1 

D!1)>SP0*SIN!B)-E!I3) 

0!2)>SP0«C0S!B)-E! 14) 

IM2-2XIMAX 
DO  52  I - 1 . I M2 
52  H!l )>0. 

H!1.I3)>1.  SH!2,l4)xl. 

CALL  R0TATE!S1.S2,B.W!1. 1 ),W!1.2).W!2.2)) 
W!2.1)«W!1,2) 

CALL  KALMAN! 2. F) 

LL«2  STM»TP  SGO  TO  110 


«»  OUTPUT 
110  CONTINUE 

KOZ«KO  SKZ>K 

IF!TP.EQ.TPS.AND.LL.EQ. 1 ) GO  TO  130 
WRITE!6,650) 

WRITE! 6. 530)  TO 
WRITE! 6, 535)  TM 
WRITE!6.700) 

CALL  UPDATE !0) 

IF!KM.EQ. 1 ) GO  TO  120 


52 


o o no 


••  OUTPUT  UNIT  LOCATION  RELATIVE  TO  OBSERVER 
WRI TE(e, 040) 

KMM-KM- I 
DO  no  K0-)  .KMM 
J2«K0*4-2  SJI«J2-1 
KK«KO+l 

DO  I 10  K«KK,KN 
l2>K*4-2  SI I•I2-I 
EX>EE( I I )-EE< Jl  ) 

EY«EE( I2)-EE( J2) 

RNO'O.  SBR0>0. 

IPtEX.EQ.O. . AND.EY.EQ.O. )00  TO  112 
RNO«SQRT(EX»EX+EY«EY) 

BR0*ATAN2(EX,EY> 

IF(BRO.LT.O. )BRO>TPtTBRO 
BRO«BRO/UB 

112  VII  -VVC  I I , I I ) -2.  •VVC  I I , Jl  I'fWI  Jl  . Jl  ) 

V22«VV( 12,  12) -2. *VV( I 2, J2)+VV( J2, J2) 

V12»VV( n , 12) -VV( 11 .J2) -VV( Jl , I2)+VV< Jl , J2) 
CALL  ELL  IPSE (VI  I ,VI2,V22,R) ,R2,AA) 

ANO'AA/UB 
CEP*FCEP(RI .R2) 

no  WRITE(6,BB0)  KO , K , RNO , BRO , Rl  , R2 , ANO , CEP 
WRITE(6.700) 

••  OUTPUT  UNIT  STATE 

120  CONTINUE 
WRITE(6,B4B) 

DO  120  K«l ,KM 

I4>K«4  SI3-I4-1  SI2-I3-1  sn>12-1 
VI  I -VV( 11,11) 

V12»VV( 11,12) 

V22>VV( 12,12) 

CALL  ELLIPSE(VI I ,VI2,V22,R1 ,R2,AA) 

ANO*AA/UB 
CEP«FCEP(R1 ,R2) 

IF(EE( 13) .EQ. 0. . AND. EE( 14) . EO. 0. )00  TO  121 
C-ATAN2(EE( 13) ,EE( 14)  ) 

IF(C. LT. 0. )C"TPI +C 
C-C/UB 
00  TO  122 

121  C-0. 

122  S>SQRT(EE(  l3)««2'fEE(  I4)*«2) 

V33«VV( 13,13) 

V34*VV( 13, 14) 

V44«VV( 14,14) 

CALL  ELLIPSE(V33,V34, V44,A1 ,A2, AA) 

AAA«AA/UB 
VEP»FCEP<A1 , A2) 

120  WRITE(6,000)  K , EE( 1 I ) , EEC  I 2) , Rl , R2, ANO, CEP, 
C,S,  A1  , A2,  AAA,VEP 
C 

130  CONTINUE 

IF(LL.EQ.2)  00  TO  0 
KO-KOZ  SK>KZ 
TPS«TP  STM-TD 

IF(TO.OT.TDS)  TDO«(TD-T0)/00. 

T0S»T0 

WRITE<6,800) 

WRITE(6,090) 

WRITE (0,600) TP, TD, KEY , KO, K . X 
00  TO  (30,40,44,00,0)  KEY 


.1  jju, 


450  FORMAT (1  HU 
600  FORMATUH  ) 

700  FORMATUH  ) 

090  FORMAT («  DATA  TP  «. 

+ » TD  KEY  KO  K XI  X2  X3  X4  XO  ») 

600  F0RMAT(6X.2F6. 1 . J3. 14, 12.  5F6. 1) 

620  F0RMAT(2F10.0. I 1 , 14, 1 2, 3X, 5F1 0, 0) 


030 

FORMAT! « 

TIME 

OF  LAST  DATA 

-*, 

F6.  1 ) 

035 

FORMAT!* 

TIME 

OF  PREDICTION 

-*, 

F6.  1) 

040 

FORMAT!* 

OBSR 

UNIT  RNQ 

BRO 

SI01 

SI02 

ANO 

CEP*) 

500 

FORMAT! 15, 

16, F7. 1 ,6F6. 1 ) 

545 

FORMAT!* 

UNIT  X 

Y 

SI01 

SIQ2 

ANO 

CEP*, 

9X. 

*CRS 

SPO 

Slot 

SI02 

ANO 

CEP«) 

305  FORMATC  Ml. 1X,6F6. 1 ,6X,6F6.  1 ) 
C 

150  CONTINUE 
END 


SUBROUTINE  UPDATE! L) 

COMMON/A/  KEY,TM,TO,KO,K,KM. IMAX.UB.TDO 

COMMON/ C/  E ( 36 ) , EE ( 36 ) , V ( 36 . 36 ) , VV ( 36 , 36 ) . A ( 36 , 36 ) 

T»(TM-T0)/60. 

IMM3*lMAX-3 
DO  10  1>1 . 1MM3,4 
A( I , 1+2)»T 
10  A( 1+1, 1+3)«T 
DO  20  1-1,1  MAX 
EE( I )-0 
DO  20  M-1 , IMAX 
20  EE( 1 )-EE( I )+A( 1 .M)«E(M) 

DO  30  1-1 , IMAX 
DO  30  J-1 . IMAX 
VV( 1 , Jl-O. 

DO  30  M> 1 , I MAX 
DO  30  N>1,IMAX 

VV( I . J)-VV< I , J)+A( 1 ,M)«V{M.N)«A{J.N) 

30  VV( J, I )-VV( 1 , J> 

IF(L.EQ.O)RETURN 
DO  70  1-1,1  MAX 
E( 1 )-EE( I ) 

DO  70  J-1, IMAX 
70  V(I,J)»VV(I.J) 

TO-TM 

RETURN 

END 
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SUBROUTINE  REPORT 

COMMON/A/  KEY.TM,TO,KO,K,KM, IMAX,UB,TDO 
COMMON/B/  D(2),W(2,2).H(2,36) 

COliMON/C/  E ( 36 ) . EE  ( 36 ) , V ( 36 . 36 ) . W ( 36,  36  ) . A ( 36 . 36 ) 

COMliON/E/  RNO,BRO,RSIO,BSIO,XE.YE 

IFITM.OT.TO)  CALL  UPDATEd) 

l2>K«4-2  SM-I2-1 

J2>KO«4-2  SJ1*J2-1 

B«BRO«UB 

IF(KEY.EQ.2)  00  TO  IS 
XObECM  )-E(J1  ) 

YDsEC 12) -E( J2) 

vii=  v(ii,ii)-2.»  v(n.ji)+  v(ji,jn 

V22>  V(l2,l2)-2.*  V(l2,J2)-f  V(J2.J2) 

V)2*  veil, 12)-  V(n,J2)-  V(J1,12)  +V<J1,J2) 

CALL  ELLIPSE! VI 1 . VI 2. V22. R1 , R2, AA) 

CAA>COS(AA) 

SAA>SIN(AA) 

XC>XD»CAA-YD«SAA 

YC«YD»CAA+XD»SAA 

BA-B-AA 

SBA<SIN(BA) 

CBA«COS(BA) 

IF(ABS(SBA) .LT. .0001 ) OO  TO  10 

IF(ABS(CBA) .LT. .0001 ) 00  TO  1 1 

TBA«SBA/CBA 

R«(R1/R2)  SRR«R«R 

Y»  < RR» YC+TBA«XC ) / < TBA« TBA+RR ) 

X«TBA*Y 

RNO^SQRTCXcX+YxY)  600  TO  12 

10  RNOaABS(YC)  600  TO  12 

11  RNO>ABS(XC) 

12  CONTINUE 
RSIO-RNO 

IB  CONTINUE 

R1<RN0<BSI0<UB 

R2<RS10 

XE-E(  J1  )-t'RNO<SIN(B) 

YE>E(  J2)-t-RNO<COS(B) 

D( 1 )>XE-E( 1 1 ) 

0(2)>YE-E( 12) 

IM2>2<IMAX 
00  25  l-l, IN2 
25  H( I )>0. 

H( 1 , 1 1 )<1 . 6H(2. I2)<1 . 

H(1.J1)--1.  6H(2.J2)>-1. 

CALL  R0TATE(R1,R2,B,W(1. 1 ) , W( 1 . 2) , W(2.  2) ) 

W(2, 1 )>W( 1 ,2) 

RETURN 

END 
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SUBROUTINE  KALMANCKA.F) 

C 

C «»  FILTER  EQUATIONS,  JAZWINSKI  P.270 
C »«  0 = V*H  » <H«V«H  + W)*«-l 
C ««  E > E 0«D 

C ««  V » ( I -0»H)»V«( I -0«H)  + 0«W*0 
C 

[ COMMON/A/  KEY.TM.TO.KO.K.KM, IMAX,UB.TDO 

F COMMON/B/  D(2),W(2.2).H(2,36) 

I COMMON/C/  E(36),EE(36),V(36,36).VV(36.36),A(36,36) 

' DIMENSION  VH(36,2),0(36.2),P(36,36),VS(36,36),U(2.2) 

C 

I4>K»4  S13=I4-1  $12=13-1  $11=12-1 
L-0 
C 

1 DO  10  1=1,1  MAX 
DO  10  J=1,2 
SUM=0. 

DO  S N« 1 , I MAX 
5 SUM=SUM-t-V(I,N)»H(J,N) 

10  VH(I,J)=SUM 
C 

DO  20  1=1,2 
DO  20  J=1,2 
SUM=W(I,  J) 

DO  15  N=1 , IMAX 
15  SUM=SUM+H(I,N)«VH(N,  J) 

20  U(I,J)=SUM 
C 

DET=U(1, 1 )»U(2,2)-U(1,2)»U(2.  1 ) 

U11=U(1,1) 

U(1 , 1 )=U(2,2)/DET 
U(1,2)=-U(1,2)/DET 
U(2,  n=U(1,2) 

U(2,2>=U1 1/DET 

IFCL.EQ. 1 .OR.TDO.EQ.O. ) 00  TO  24 
C 

R=0. 

DO  21  M=1,2 
DO  21  N=1,2 

21  R=R-t-0(M)«U(M,N)«0(N) 

F=EXP(-R/2.)  $FF=1.-F 

IF(KA.EQ.I)  RETURN 

C 

00  TO  (22,22,22,23)  KEY 

22  V(I1.I1)=V(I1,I1)  * 0(1)«D(1)«FF 
V( 12, I2)=V( 12, 12)  + D(2)»D(2)«FF 
V(I1,  I2)«V(I1,  12)  -I-  D(1)«D(2)»FF 
V(I2.I1)=V(11,I2) 

D3=D(1 )/TDO 
D4=0(2)/T0e 

V( 13, I3)=V( 13,  13)  + 03»03»FF 
V(14, I4)=V(I4, 14)  * D4«04«FF 
V( 13, I4)=V( 13, 14)  * D3«D4«FF 
V( 14, I3)=V( 13,  14) 

V(l  1,  13)=V(  1 1,  13)  D(1)«03«FF 

V( 1 1 , I4)=V( 1 1 , 14)  * D(1)*D4«FF 
V(  12,  I3)=V(  12,  13)  D(2)«D3*FF 

V(I2, I4)=V(I2, 14)  * D(2)«04«FF 
V(13,I1)=V(I1.I3) 

V(I4,I1)=V(I1,I4) 

V( 13, I2)=V( 12, 13) 

V(14,I2)=V(I2,I4) 

L=1  $00  TO  1 
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c 


c 


c 


c 


23  V(  13,  I3)>V(  13.  13)  D(1)«D(1)«FF 

V(  14,  I4)«V(  14.  14)  -f  0(2)»D(2)«FF 
V(  13.  I4)«V(  13,  14)  0(1)S0(2)«FF 

V( 14, I3)«V( 13. 14) 

L>1  SOe  TO  1 

24  DO  30  I « 1 , I MAX 
DO  30  J-1 ,2 
SUM>0. 

DO  25  N«1,2 

25  SUM-SUM-)-VH(l,N)«U(N,  J) 

30  0( I , J)>SUM 

DO  40  I -1 , IMAX 
DO  40  , IMAX 

SUM^O. 

IF( I .EQ. J)SUM»1  . 

DO  35  N-1,2 

35  SUM>SUM-0( I ,N)»H(N, J) 

40  P( I , J)3SUM 

DO  55  1 = 1,1  MAX 
DO  55  J=1 , IMAX 
SUM=0. 

DO  45  N=1, IMAX 
DO  45  M= 1 , 1 MAX 

45  SUM»SUM+P( 1 ,N)*V(N.M)»P(J.M) 

DO  50  N=1,2 
DO  50  M- 1 , 2 

50  SUM»SUM+0( 1 ,N)»W(N,M)*D( J.M) 

55  VS( 1 , J)«VS( J, 1 )=SUM 

DO  65  1=1, IMAX 
DO  60  N-1,2 

60  E(l)>E(l)-«-0(l,N)«0(N) 

DO  65  J«1 , IMAX 
65  V(1,J)»V(J.1)»VS(1,J) 

RETURN 

END 
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SUBROUTINE  ELLIPSE! VI 1 , VI 2. V22. R1 . R2. A) 
IP(V12.EQ.O. ) 00  TO  10 
A>.5«ATAN2(2.sV12.V22-V1 1 ) 

SINA>SIN(A) 

COSA>COS(A) 

SSIN>SINA*SINA 

SCOS>SINA«COSA 

CCOS«COSA«COSA 

RR1  >CC0S«V1 1 -2.  «SC0S«V12-fSSIN«V22 
RR2>SSIN«V1 1-I-2.  «SC0S«V12-t-CC0S»V22 
R1«SQRT(RR1 ) 

R2>SQRT(RR2) 

RETURN 

10  R1<SQRT(V1 1 ) 

R2>SQRT(V22) 

A*0. 

IF(R2.0E.R1)  RETURN 

RbRI  SR1sR2  SR23R  SA>3 . 141 S926S359/2 . 

RETURN 

END 


SUBROUTINE  ROTATEIRI . R2, A. VI 1 , V12.V22) 
SINA>S1N(A) 

COSA^COSCA) 

SSIN>SINA«SINA 

SCOS>SINA«COSA 

CC0S3C0SA«C0SA 

RR1>R1«R1 

RR2>R2»R2 

V1 1 ■CC0S«RR1  -l-SSI  N«RR2 

V12«SC0S«(RR2-RR1 ) 

V22=SSIN»RR1-«-CC0S*RR2 

RETURN 

END 


FUNCTION  FCEP(R1,R2) 

CEP- . S62«R2-t- . 6 1 0-RI 

IFIRI  .LT.  .3«R2)  CEP-R2*  (.  678-t- . 79444*  ( R1 /R2>  ««2) 

FCEP-CEP 

RETURN 

END 
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